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We apply a Color Glass Condensate+Non-Relativistic QCD (CGC+NRQCD) framework to com¬ 
pute J/’if production in deuteron-nucleus collisions at RHIC and proton-nucleus collisions at the 
LHC. Our results match smoothly at high to a next-to-leading order perturbative QCD -|- 
NRQCD computation. Excellent agreement is obtained for px spectra at RHIC and LHC for cen¬ 
tral and forward rapidities, as well as for the normalized ratio RpA of these results to spectra in 
proton-proton collisions. In particular, we observe that the RpA data is strongly bounded by our 
computations of the same for each of the individual NRQCD channels; this result provides strong 
evidence that our description is robust against uncertainties in initial conditions and hadronization 
mechanisms. 

PACS numbers: 11.80.La, 12.38.Bx, 14.40.Pq 


The copious production of heavy quarkonium states at 
high energy colliders has inaugurated a new era of preci¬ 
sion studies of such states [I| . In proton-proton collisions 
(p-bp), next-to-leading order (NLO) perturbative stud¬ 
ies are available within the Non-Relativistic QCD 
(NRQCD) factorization framework These computa¬ 
tions can be further improved by employing QCD fac¬ 
torization [3, 0 to resum large logarithms ln(px/AL) in 
the ratio of the transverse momentum px to the quark 
mass M. A comparison of these studies with collider 
data therefore provides key insight into the formation 
and hadronization of heavy quark-antiquark pair {QQ- 
pair) states in QCD. 

In proton-nucleus (p+A) collisions, additional features 
of QQ-pair production and hadronization can be tested. 
These include many-body QCD effects such as multiple 
scattering and shadowing of gluon distributions in nuclei, 
as well as the radiative energy loss induced in the scatter¬ 
ing of the QQ-pair off the colored medium. Besides these 
insights into many-body QCD dynamics, p-|-A collisions 
also provide an important benchmark for understanding 
the interactions of heavy quarks in the hot and dense 
medium created in heavy ion collisions. 

For small gluon momentum fractions x, their distri¬ 
butions saturate with a dynamically generated satura¬ 
tion scale Qs{x) [im. This regime is accessed when 
P± ^ Qs- The Color Glass Condensate (CGC) effec¬ 
tive theory 0111 provides a quantitative framework 
to study many-body QCD effects in high energy scat¬ 
tering processes when Qs(x) » Aqcd; where Aqcd is 
the fundamental QCD scale. In this limit, multiple scat¬ 
tering contributions can be absorbed into light like Wil¬ 
son line correlators, which govern the shadowing and px 
broadening of QQ-pair distributions at small x. Energy 
evolution of these correlations at small x is described 
by the Balitsky-JIMWLK hierarchy of renormalization 
group equations [iM3. Energy loss contributions, in¬ 
cluded in some models in the literature [l^ , are formally 


NLO in the CGC framework (T^ . 

Expressions for QQ-pair production in p-l-A collisions 
in the CGC framework were derived previousl y in [ 13 - 
1^ as well as in related dipole approaches For 

Px >> Qs, the results can be matched to those derived 
in perturbative QCD frameworks [1^. In [s^, the ma¬ 
trix elements in the CGC framework were combined with 
the Color Evaporation hadronization Model (CEM) to 
compute J/if production in proton-proton and proton 
(deuteron)-nucleus collisions at the LHC (RHIC)^. The 
quantity 


RpA — 


d(JpA 
A X dap 


( 1 ) 


the ratio of the cross-sections in p-bA collisions to p+p 
collisions, normalized by the atomic numbe r A, was found 
to be suppressed relative to the data [3l|, ■ Very re¬ 

cently, the authors of [s^ argued that better agreement 
of the CGC-bCEM model with the RpA data is obtained 
if nuclear effects were treated differently. Here we shall 
apply NRQCD to describe the hadronization of QQ-pair 
and compute J/'tf production in a CGC-bNRQCD frame¬ 
work [34| . In addition to providing a more systematic 
power counting, NRQCD allows one to smoothly match 
the CGC computations to successful NLO NRQCD com¬ 
putations for Px >> Qs- This strategy was previously 
applied to successfully describe p-bp collisions at RHIC 
and LHC [H. 

Eor c omp leteness, we outline the CGC-bNRQCD for¬ 
malism [^[ 33 . In NRQCD factorization, the production 
cross section of a quarkonium H in the forward region of 


^ For simplicity, we will generically call both sorts of collisions p+A 
collisions in the rest of the paper. 
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a p+A collision is expressed as @ 

= ( 2 ) 

K 

where k = denotes the quantum numbers of the 

intermediate QQ-pair in the standard spectroscopic no¬ 
tation for angular momentum. The superscript c denotes 
the color state of the pair, which can be either color sin¬ 
glet (CS) with c = 1 or color octet (CO) with c = 8. 
For J/'i/j production that will be studied here, the most 
important intermediate states are and 

In Eq. ([2]), (O^) are non-perturbative universal 
long distance matrix elements (LDMEs), which can be 
extracted from data, and dd” are short-distance coeffi¬ 
cients (SDCs) for the production of a QQ-pair, computed 
in perturbative QCD. 

To calculate the SDCs in Eq. ([2]), we apply the CGC 
effective field theory which results in the expres¬ 

sions [35l|. 

CS cis{'n'R\) f ^p,yp (fcm) 

d^p±dy (27r)9(Af2 _ X) J 

X My ik±)J\fY{k'j_)Jd'Y{p± — feij_ — fci. — 
for the color-singlet ^ channel, and 

CO otsiirR'^) f ^p,yp{ki±) 

d^PYdy (27r)7(7V2 - 1) J k‘fj_ 

ki±,k± 

X A/V(fcj.)A/V(pj. — ki± — fej_) Tg, 

for the color-octet channels. Here <Pp,j/p is the uninte¬ 
grated gluon distribution inside the proton, which can 
be expressed as 

_ AT 1.2 ^ 

<Pp,yp{ki±) = TrRl-^^My^{ki±). (5) 

The functions Qi (Tg) are calculated perturbatively-the 
expressions are available in [H (dl). M {M^) are the 
momentum-space dipole forward scattering amplitudes 
with Wilson lines in the fundamental (adjoint) represen¬ 
tation, and TrRp {'kR\) is the effective transverse area 
of the dilute proton (dense nucleus). These formulas can 
be used to compute quarkonium production in p-l-A colli¬ 
sions. By replacing “A’s by p’s”, they can also be used to 
compute quarkonium production in p+p collisions [s^. 
For deuteron-gold (d-|-Au) collisions at RHIC, since gluon 
shadowing effects are weak for deuteron side, we assume 

= 2 (/3p^yp(feij_). 

Before we confront our framework to data on p-l-A col¬ 
lisions, there are a number of parameters that have to 
be fixed. Nearly all the parameters are identical to those 
previously determined in our study [s^ of p-|-p collisions. 
The charm quark mass is set to be m = 1.5 GeV, approx¬ 
imately one half the J/ijr mass. The CO LDMEs were ex¬ 
tracted in the NLO collinear factorized NRQCD formal¬ 
ism Q by fitting Tevatron high p_L prompt J/ip produc¬ 
tion data; one obtains {O'^M= 1.16/(2iVc) GeV^, 


^ 0 j/V-(i 5 '[ 8 l)) = 0.089 ± 0.0098 GeV^ {O'’M= 
0.0030 ± 0.0012 GeV^ and {O-^M = 0.0056 ± 
0.0021 GeV^. We emphasize, as previously, that the high 
sensitivity of short distance cross-sections to quark mass 
can be mitigated by the mass dependence of the LDMEs. 
Note that the uncertainties of these CO LDMEs include 
only uncorrelated statistic errors, but not correlated er¬ 
rors Q. Further, as in [s^, M and M"^ are obtained by 
solving the running coupling Balitsky-Kovchegov (rcBK) 
equation [TsI . in momentum space with McLerran- 
Venugopalan (MV) initial conditions [ll|> B Ibe 
dipole amplitude at the initial rapidity scale Yq = 
ln(l/a:o) (with xq = 0.01) for small x evolution. In the 
case of p-l-p collisions, all the parameters in the rcBK 
evolution are fixed from fits to the HERA DIS data [s^ ■ 
In [^, we devised a matching scheme that allowed us 
to interpolate between the proton’s collinearly factorized 
gluon distribution at large x with the unintegrated dis¬ 
tribution in Eq. ©• This allowed us to fix the remaining 
free parameter, the effective gluon radius of the proton 
Rp = 0.48 fm. 

Turning to p-|-A collisions, there are two additional pa¬ 
rameters in our framework, the initial saturation scale 
Qso,a in the nucleus and the effective transverse ra¬ 
dius Ra- The latter is not the charge radius of the 
nucleus, but parametrizes the overall non-perturbative 
cross-section of relevance to quarkonium production. A 
more detailed treatment would take into account the im¬ 
pact parameter dependence of the unintegrated distri¬ 
butions, and model the inelastic proton-nucleus cross- 
section as in B, . We will return to this point shortly. 
In general, we can express the initial saturation scale in 
the nucleus as = Mx Qgg.pi where V is a number to 

be determined and Q^ p is the initial saturation scale in 
proton, fixed by the fit to HERA DIS data [s^ ■ Good fits 
to extant electron-nucleus (e-l-A) DIS data were obtained 
in for rcBK evolution with the following initial condi¬ 
tions: i) MV model with anomalous dimension 7 = 1.13, 
ii) MV model with anomalous dimension 7 = 1 . For the 
initial conditions i), one obtains a good fit to e-|-A data 
for N 3, while for initial conditions ii), N ss 1.5. In 
this paper, rcBK evolution for nuclei was performed for 
initial conditions ii). To avoid fine tuning, we will choose 
N = 2 for the results presented in this paper^. 

Similar to Rp for the proton, the effective radius 
Ra providing the non-perturbative normalization of the 
cross-section here can be different from the transverse 
charge radius of the nucleus because we have a specific 
heavy particle produced in the final state. Fortunately, 
there is a physical condition which we can use to con- 


^ The quality of fit for Af = 2 is marginally better than that for 
N = 3 but significantly better than those for Af = 1 or higher 
values of integer N. For the IP-sat model, for median impact 
parameters in e-|-A DIS, = 6 in contrast to A ss 2 for b = 

0 [ 43 . 
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strain it. When pj_ is much larger than the saturation 
scale involved, the gluon distribution becomes dilute and 
the nuclear suppression effect should be negligible. Thus 
RpA must approach 1 for high p±_ » Qso.a- Using 
Eqs. (P)-® one can derive (the argument is presented 
in Appendix), the expression. 


R\ QVo.a 

AH qX 


( 6 ) 


We will see later that Eq. ® indeed guarantees RpA 1 
at highp_L limit, within a few percent. By choosing 7=1 
and N = 2, we obtain Ra = \fAj2Rp^ which equals to 
4.9 fm for Pb and 4.8 fm for Au^. 

Because Eqs. are computed only at LO in 

the CGC power counting, the CGC+NRQCD framework 
cannot be extended to describe high p± p+p and p+A 
data, one might challenge that using Eq. ® to deter¬ 
mine Ra is not especially meaningful. We emphasize 
however that this condition must be satisfied for the 
GGG-I-NRQGD framework to be self-consistent at each 
order in the perturbative expansion. The p± at which 
Eq. ® is saturated may differ. At NLO, the above proce¬ 
dure should be redone to determine a new self-consistent 
condition. 


Vs=5.02TeV for LHC pPb 
Vs=0.2TeV for RHIC dAu 


mm CGC+NRQCD 
» NLONRQCD 
, LHCb data for 1.5<y<4.0 

* LHCb data for 1.5<y<2.0 

, 0.3xLHCh data for 2.0<y<2.5 
, O.lxLHCb data for 2.5<y<3.0 

* 0.03xLHCb data for 3.0<y<3.5 
. O.OlxLHCb data for 3.5<y<4.0 
^ 0.003XRHIC data for 1.2<y<2.2 
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FIG. 1: p± spectrum of J/i/> production in p+Pb collisions 
at 5.02 TeV and d+Au collisions at 0.2 TeV. NLO NRQCD 
results are taken from Ref. [4^ . The experimental data are 
taken from Refs. [^. [^. 


To better present the p-|-A results, we define a cross 
section per nucleon-nucleon collision, daNN = 


® Interestingly, the ratio Ra/Rp ~ 10 here is close to the ratio 
of radii extracted from estimates of the inelastic p+A and p+p 
cross-sections at both LHC and RHIC. 



y 


FIG. 2: Rapidity distribution of J/'ip production in p+Pb 
collisions at 5.02 TeV and d+Au collisions at 0.2 TeV. The 
experimental data are taken from Refs, 


Fig. [T] displays the p^ spectrum of J/t/> production in 
p-|-Pb collisions at 5.02 TeV and d-|-Au collisions at 
0.2 TeV. The bands of our GGG-I-NRQCD results es¬ 
timate uncorrelated errors of LDMEs and an additional 
global 30% uncertainty to account for correlated errors 
of LDMEs, errors from treatment of feed down, veloc¬ 
ity corrections and radiative corrections. We find that 
the contribution of the GS channel is about 15 — 20% at 
small p± and decreases as pj_ becomes larger. The NLO 
NRQGD predictions are taken from [d^, where the PDF 
shadowing model EPS09 [dj was employed to estimate 
the (small) nuclear shadowing effects at large p±. For all 
rapidity bins available, the GGG-I-NRQGD curves match 
on to the NLO NRQGD ones smoothly, providing a good 
description of all experimental data. Interestingly, one 
finds that the GGG+NRQGD curves overshoot the data 
at smaller values of p± at RHIG relative to the LHG data. 
This may be anticipated because, for a given p±, small x 
logs are less important at lower energies. However, a full 
NLO computation in this framework is needed to under¬ 
stand better the matching in p± of the two formalisms. 

The rapidity distribution of J/^/> production in p+Pb 
collisions at 5.02 TeV and d-fAu collisions at 0.2 TeV 
is shown in Fig. [5J where the the bands are generated 
similarly to those in Fig. [TJ Since these data are inte¬ 
grated over p±, the low p± region dominates and the 
GGG+NRQGD formalism at LO should apply. Both 
LHG data and forward RHIG data are well covered by 
our uncertainty band; the central value for mid-rapidity 
RHIG data however is slightly below the band. For this 
data point, our theory curves should have a larger sys¬ 
tematic uncertainty because our framework is most re¬ 
liable for dilute-dense collisions corresponding to high 
energies and forward rapidities. The key observation 
though is that both the relative shapes as well as the 
absolute magnitudes of the curves are well captured in 
the GGG+NRQGD formalism. The quality of the fits to 
the pj_ and rapidity spectra in Figs. [1] and [5] are similar 
to those in p+p collisions [s^. Thus we should be able 
to describe the RpA ratio, which we shall now discuss. 

A key point is that the the large uncertainties for 
LDMEs, feed down contributions and velocity correc¬ 
tions, largely cancel in the ratio of each NRQGD channel 


















4 





FIG. 3: RpA as a function of p± (upper) and rapidity (lower) 
at LHC. The experimental data are taken from Refs. [^. [^. 

US- 


FIG. 4: RpA as a function of p± (upper) and rapidity (lower) 
at RHIC. The experimental data are taken from Refs. [^.IdSl. 


contributing to J/'^ production. The band spanned by 
different channels should be able to bracket the RpA value 
for J/ijj production. With this method, the bounded 
value of RpA extracted for J/?/; production is independent 
of the LDMEs and their statistical uncertainties. This 
is especially noteworthy since independent extractions of 
the LDMEs from present data are not feasible; their mag¬ 
nitudes, especially between the various CO channels, can 
vary significantly. Finally, since the CEM is a special 
case of NRQCD with the choice of certain LDMEs [d^ , 
our calculation of RpA will also cover the range of CEM 
predictions. In this sense, the range of theoretical esti¬ 
mates of RpA for J/il) production are independent of the 
J/'0 hadronization model and are directly sensitive to the 
short distance physics. 

We will employ here the principal channels for J/ip 
production given by NRQCD power counting-these cor¬ 
respond to the and channels. Our 

results for RpA as a function of p± and rapidity, com¬ 
pared to data from the LHC and RHIC, respectively, are 
presented in Figs. |3] and |4l where a 5% systematical error 
is assumed for each channel to account for the approx¬ 
imation in Eq. The RpA of all CO channels ap¬ 

proaches 1 at high pj_ , confirming that condition Eq. 
indeed is satisfied by the full theoretical calculation. On 
the contrary, RpA of the CS channel increases to 
be larger than 1 at high p±. Since forming a color sin¬ 
glet requires two gluons from the target, the additional 
gluon exchange from the nucleus, at high pj _, is enhanced 
relative to that from a proton (by an amount that is pro¬ 
portional asymptotically to the ratio of their saturation 
scales at the rapidity of interest). Nevertheless, as we 
find the contribution of the CS channel is small relative 
to the CO terms in both p-l-p and p-l-A collisions, it does 
not affect our estimate of RpA- Thus the band represent¬ 


ing the RpA spanned by the CO channels corresponds to 
our result for RpA oi J/ip production. 

The p± and rapidity RpA data from both RHIC and 
LHC lie within our uncertainty bands. At the LHC, the 
35(8] 

state lies closest to the central values of the data, 
while at RHIC, the and ^Pj®^ channels are closest 
to the data. Our results suggest that the RpA data, in 
a future global analysis within the CGC/NLO-I-NRQCD 
framework, can help constrain the LDMEs more strin¬ 
gently, thereby providing a further test of NRQCD. 

To summarize, we have shown here that J/ip spectra 
in p-l-A collisions both at RHIC and the LHC are well 
described by our CGC-I-NRQCD computations. The two 
free non-perturbative parameters are related by Eq. ®; 
further, the value of the initial nuclear saturation scale 
Qso,a is consistent with the values that best describe 
fixed target e-l-A DIS data. The fact that the RpA p± 
data lie within the bands spanned by our computations 
for the different color octet channels is a strong evidence 
for the robustness of our framework since these curves are 
insensitive to details of how heavy quark pairs hadronize 
to form the J/ip- The results in this paper, when com¬ 
bined with those in (35j| . provide the first comprehensive 
description of J/ip production in both p-l-p and p-l-A col¬ 
lisions at collider energies. 

Several outstanding questions remain. Firstly, the 
NLO CGC computation needs to be performed to confirm 
that the framework established is solid. Secondly, other 
quarkonium states remain to be studied; these come with 
unique challenges. For instance, for T production, Su- 
dakov type double lo gs in M/P± are important and need 
to be resummed [48l - l50l |. A systematic computation of 
ip(2S) production in p-l-A collisions, may require that we 
include relativistic contributions in the computation of 
the heavy quark matrix elements. All these questions 
can be explored in the framework discussed here. 
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Appendix A: Derivation of Eq. ([B]) 


decreases as inverse powers of fc_L, the dominant 
contribution for Eq. comes from two regions, either 
fcj_ is small or p± — k± is small. It is clear that Eg 
must be a symmetric function under the transformation 
k± —>■ p± — k±, thus, in either case, we can set k± to be 
zero in Eg and it becomes independent of k±. Then we 
can perform the k± integration in Eq. dsu, which gives 
Therefore, the RpA defined in Eq. m behaves 
as 


Let us derive the corresponding relation from this con¬ 
dition. As CS contribution is negligible at very high p± 
regime [s^, we will only consider the CO contribution in 
Eq. dH). When p± is large, at least one of ki±, k± and 
p± — ki± — k± needs to be large. As we are considering 
a dilute-dense collision, the contribution from the region 
where fci_L is large should be less important, which im¬ 
plies we can take the collinear limit for the proton side 
and give [s^ 


da^ 

d?pYdy 


CO 


as{TTR\) 


4(27r)3(iV2 - 1) 


^pfp/gi^pjQ ) 


X y’w(fex)W(px-/cj.)f^, 


(Al) 


where fp/g is the gluon collinear PDF and E g ar e collinear 
limit of Eg which have been calculated in |^ . Because 




Ra J^Ya^Pi^ 


(A2) 


Note that, although we have used the collinear ap¬ 
proximation to derive the above relation, the relation 
holds much better than the collinear approximation itself, 
which is caused by the cancellation between the contri¬ 
butions to the numerator and denominator of RpA from 
large fcij_ region. It is known that MV model with rcBK 
equation gives A/"^ (pj_) oc QYa high pj_ limit [l^, we 


therefore have Y^a^Y^] 
J^d{p±) 


S; s ,p 


—where at the last 


VsO,p 


step we assume that the Y is not significantly larger than 
Yq and thus the ratio of saturation scales is not changed 
too much by evolution. Following these steps, we obtain 
the condition in Eq. 


[1] N. Brambilla, S. Eidelman, B. Heltsley, R. Vogt, 

G. Bodwin, et ai, Heavy quarkonium: progress, puzzles, 
and opportunities, Eur.Phys.J. C71 (2011) 1534 
arXiv:1010.5827 InSPIRE . 

[2] M. Butenschoen and B. A. Kniehl, J/’ijj polarization at 
Tevatron and LHC: Nonrelativistic-QCD factorization 
at the crossroads, Phys.Rev.Lett. 108 (2012) 172002 
arXiv:1201.1872 InSPIRE . 

[3] K.-T. Chao, Y.-Q.^ Ma, H.-S.^ Shao, K. Wang, and Y.-J. 
Zhang, J/%p Polarization at Hadron Colliders in 

Nonrelativistic QCD, Phys.Rev.Lett. 108 (2012) 242004 
arXiv:1201.2675 InSPIRE . 

[4] B. Gong, L.-R Wan, J.-X. Wang, and H.-F. Zhang, 

Polarization for Prompt J/ip, production at the 

Tevatron and LHC, Phys.Rev.Lett. 110 (2013) 042002 
arXiv:1205.6682 InSPIRE . 

[5] H. Shao, H. Han, y! Ma, C. Meng, Y. Zhang, et ai, 
Yields and polarizations of prompt J/tf and ‘4i{2S) 
production in hadronic collisions, arXiv: 1411.3300 

InSPIRE . 

[6] G. T. Bodwin, E. Braaten, and G. P. Lepage, Rigorous 
QCD analysis of inclusive annihilation and production 
of heavy quarkonium, Phys.Rev. D51 (1995) 1125-1171 
hep-ph/9407339 InSPIRE . 

[7] Z.-B. Kang, Y.-Q. Ma, J.-W. Qiu, and G. Sterman, 
Heavy quarkonium production at collider energies: 
Factorization and Evolution, 

Phys.Rev. D90 (2014) 034006 arXiv: 1401.0923 


InSPIRE . 

[8] Z.-B. Kang, Y.-Q. Ma, J.-W. Qiu, and G. Sterman, 
Heavy Quarkonium Production at Collider Energies: 
Partonic Cross Section and Polarization, 

Phys.Rev. D91 (2015) 014030 arXiv: 1411.2456 
InSPIRE . 

[9] L. Gribov, E. Levin, and M. Ryskin, Semihard Processes 
in QCD, Phys.Rept. 100 (1983) 1-150 InSPIRE . 

[10] A. H. Mueller and J. Qiu, Cluon Recombination and 
Shadowing at Small Values of x, 

Nucl.Phys. B268 (1986) 427 InSPIRE . 

[11] L. D. McLerran and R. Venugopalan, Computing quark 
and gluon distribution functions for very large nuclei, 
Phys.Rev. D49 (1994) 2233-2241 hep-ph/9309289 
InSPIRE . 

[12] L. D. McLerran and R. Venugopalan, Cluon distribution 
functions for very large nuclei at small transverse 
momentum, Phys. Rev. D49 (1994) 3352-3355 
hep-ph/9311205 InSPIRE . 

[13] E. lancu and R. Venugopalan, The Color glass 
condensate and high-energy scattering in QCD, 
hep-ph/0303204 InSPIRE . 

[14] F. Gelis, E. lancu, J. Jalilian-Marian, and 

R. Venugopalan, The Color Class Condensate, 

Ann. Rev. Nucl. Part. Sci. 60 (2010) 463-489 
arXiv: 1002.0333 InSPIRE. 

[15] 1. Balitsky, Operator expansion for high-energy 
scattering, Nucl.Phys. B463 (1996) 99-160 









6 


hep-ph/9509348 InSPIRE . 

[16] E. lancu, A. Leonidov, and L. D. McLerran, Nonlinear 
gluon evolution in the color glass condensate. 1., 
Nucl.Phys. A692 (2001) 583-645 hep-ph/0011241 
InSPIRE . 

[17] J. Jalilian-Marian, A. Kovner, and H. Weigert, The 
Wilson renormalization group for low x physics: Gluon 
evolution at finite parton density, 

Phys.Rev. D59 (1998) 014015 hep-ph/9709432 
InSPIRE . 

[18] F. Arleo and S. Peigne, Heavy-quarkonium suppression 
in p-A collisions from parton energy loss in cold QCD 
matter, J77EP 1303 (2013) 122 arXiv: 1212.0434 
InSPIRE . 

[19] T. Lion and A. Mueller, Parton energy loss in high 
energy hard forward processes in proton-nucleus 
collisions, Phys.Rev. D89 (2014) 074026 

arXiv: 1402.1647 InSPIRE. 

[20] J. P. Blaizot, F. Cells, and R. Venugopalan, High-energy 
pA collisions in the color glass condensate approach. 2. 
Quark production, Nucl.Phys. A743 (2004) 57-91 
hep-ph/0402257 InSPIRE . 

[21] H. Fuji!, F. Cells, and R. Venugopalan, Quantitative 
study of the violation of k±-factorization in 
hadroproduction of quarks at collider energies, 

Phys.Rev.Lett. 95 (2005) 162002 hep-ph/0504047 
InSPIRE . 

[22] H. Fuji!, F. Cells, and R. Venugopalan, Quark pair 
production in high energy pA collisions: General 
features, Nucl.Phys. A780 (2006) 146-174 
hep-ph/0603099 InSPIRE . 

[23] D. Kharzeev and K. Tuchin, Signatures of the color 
glass condensate in J/%p production off nuclear targets, 
Nucl.Phys. A770 (2006) 40-56 hep-ph/0510358 
InSPIRE . 

[24] D. Kharzeev, E. Levin, M. Nardi, and K. Tuchin, Gluon 
saturation effects on J/if production in heavy ion 
collisions, Phys.Rev.Lett. 102 (2009) 152301 

arXiv: 0808.2954 InSPIRE . 

[25] F. Dominguez, D. Kharzeev, E. Levin, A. Mueller, and 
K. Tuchin, Gluon saturation effects on the color singlet 
J /if) production in high energy dA and AA collisions, 
Phys.Lett. B710 (2012) 182-187 arXiv: 1109.1250 
InSPIRE . 

[26] E. Akcakaya, A. Schafer, and J. Zhou, Azimuthal 
asymmetries for quark pair production in pA collisions, 
Phys.Rev. D87 (2013) 054010 arXiv: 1208.4965 
InSPIRE . 

[27] B. Kopeliovich and A. Tarasov, Gluon shadowing and 
heavy flavor production off nuclei, 

Nucl.Phys. A710 (2002) 180-217 hep-ph/0205151 
InSPIRE . 

[28] L. Motyka and M. Sadzikowski, On relevance of triple 
gluon fusion in J/if hadroproduction, 

arXiv: 1501.04915 InSPIRE. 

[29] F. Cells and R. Venugopalan, Large mass qq production 
from the color glass condensate, 

Phys.Rev. D69 (2004) 014019 hep-ph/0310090 
InSPIRE . 

[30] H. Fuji! and K. Watanabe, Heavy quark pair production 
in high energy pA collisions: Quarkonium, 

Nucl.Phys. A915 (2013) 1-23 arXiv: 1304.2221 
InSPIRE . 

[31] ALICE Collaboration , B. B. Abelev et al., J /if 


production and nuclear effects in p-Pb collisions at 
y/S/f/f = 5.02 TeV, JHEP 1402 (2014) 073 
arXiv: 1308.6726 InSPIRE. 

[32] LHCb , R. Aaij et al., Study of J/if production and 
cold nuclear matter effects in pPb collisions at 
,/s//// =5 TeV, JHEP 1402 (2014) 072 

arXiv: 1308.6729 InSPIRE. 

[33] B. Ducloue, T. Lappi, and H. Mantysaari, Forward J/if 
production in proton-nucleus collisions at high energy, 
arXiv:1503.02789 InSPIRE . 

[34] Z.-B. Kang, Y.-Q. Ma, and R. Venugopalan, 

Quarkonium production in high energy proton-nucleus 
collisions: CGC meets NRQCD, JHEP 1401 (2014) 056 
arXiv: 1309.7337 InSPIRE. 

[35] Y.-Q. Ma and R. Venugopalan, Gomprehensive 
Description of J/if Production in Proton-Proton 
Gollisions at Collider Energies, 

Phys.Rev.Lett. 113 (2014) 192301 arXiv: 1408.4075 
InSPIRE . 

[36] Y. V. Kovchegov, Small-x F 2 structure function of a 
nucleus including multiple pomeron exchanges, 

Phys.Rev. D60 (1999) 034008 hep-ph/9901281 
InSPIRE . 

[37] J. L. Albacete, A. Dumitru, H. Fujii, and Y. Nara, CGC 
predictions for p-hPb collisions at the LHC, 

Nucl.Phys. A897 (2013) 1-27 arXiv: 1209.2001 
InSPIRE . 

[38] M. L. Miller, K. Reygers, S. J. Sanders, and 

P. Steinberg, Glauber modeling in high energy nuclear 
collisions, Ann.Rev.Nucl.Part.Sci. 57 (2007) 205-243 
nucI-ex/0701025 InSPIRE. 

[39] D. C. d’Enterria, Hard scattering cross-sections at LHC 
in the Glauber approach: From pp to pA and AA 
collisions, nucl-ex/0302016 InSPIRE . 

[40] K. Dusling, F. Celis, T. Lappi, and R. Venugopalan, 
Long range two-particle rapidity correlations in A-hA 
collisions from high energy QCD evolution, 

Nucl.Phys. A836 (2010) 159-182 arXiv:0911.2720 
InSPIRE . 

[41] H. Kowalski, T. Lappi, and R. Venugopalan, Nuclear 
enhancement of universal dynamics of high parton 
densities, Phys.Rev.Lett. 100 (2008) 022303 

arXiv:0705.3047 . 

[42] H.-F. Zhang et al. In preparation . 

[43] PHENIX Collaboration , A. Adare et al, 

Transverse-Momentum Dependence of the J/if Nuclear 
Modification in d-\-Au Collisions at /snn = 200 GeV, 
Phys.Rev. C87 (2013) 034904 arXiv: 1204.0777 
InSPIRE . 

[44] k. Eskola, H. Paukkunen, and C. Salgado, EPS09: A 
New Generation of NLO and LO Nuclear Parton 
Distribution Functions, JHEP 0904 (2009) 065 
arXiv: 0902.4154 InSPIRE. 

[45] PHENIX Collaboration , A. Adare et al.. Cold 
Nuclear Matter Effects on J/if Yields as a Function of 
Rapidity and Nuclear Geometry in Deuteron-Gold 
Collisions at /snn = 200 GeV, 

Phys.Rev.Lett. 107 (2011) 142301 arXiv: 1010.1246 
InSPIRE . 

[46] C. T. Bodwin, E. Braaten, and J. Lee, Comparison of 
the color-evaporation model and the NRQCD 
factorization approach in charmonium production, 
Phys.Rev. D72 (2005) 014004 hep-ph/0504014 
InSPIRE . 






7 


[47] ALICE , J. Adam et al, Rapidity and 
transverse-momentum dependence of the inclusive J/ip 
nuclear modification factor in p-Pb collisions at 

TeV, arXiv: 1503.07179 InSPIRE . 

[48] E. L. Berger, J. Qiu, and Y. Wang, Transverse 
momentum distribution of T production in hadronic 
collisions, Phys.Rev. D71 (2005) 034007 
hep-ph/0404158 InSPIRE . 

[49] P. Sun, C.-P. Yuan, and F. Yuan, Heavy Quarkonium 


Production at Low Pt in NRQCD with Soft Gluon 
Resummation, Phys.Rev. D88 (2013) 054008 
arXiv: 1210.3432 InSPIRE. 

[50] J.-W. Qiu, P. Sun, B.-W. Xiao, and F. Yuan, Universal 
Suppression of Heavy Quarkonium Production in pA 
Collisions at Low Transverse Momentum, 

Phys.Rev. D89 (2014) 034007 arXiv: 1310.2230 
InSPIRE . 


